On 
On 



23 



c3 



c3 



Estimating the J function without 
edge correction* 

A. Baddeley 1 , M. Kerschei 4 , K. Schladitz+, B. T. Scott 



Department of Mathematics and Statistics, University of Western Australia, 
O ■ Nedlands WA 6907, Australia 

+ Institut fur Techno- und Wirtschaftsmathematik, 
Erwin-Schrodinger-StraBe, D-67663 Kaiserslautern, Germany 

C/3 ' * Sektion Physik, Universitat Miinchen, 

^ ■ Theresienstr. 37/111, 80333 Miinchen, Germany 

The interaction between points in a spatial point process can be measured by its empty 
space function F, its nearest-neighbour distance distribution function G, and by combi- 
nations such as the J-function J = (1 — G)/(l — F). The estimation of these functions is 
hampered by edge effects: the uncorrected, empirical distributions of distances observed 
in a bounded sampling window W give severely biased estimates of F and G. How- 
ever, in this paper we show that the corresponding uncorrected estimator of the function 
J = (1 — G)/(l — F) is approximately unbiased for the Poisson case, and is useful as a 
summary statistic. Specifically consider the estimate Jw of J computed from uncorrected 



estimates of F and G. The function Jyy(r), estimated by J\y, possesses similar properties 
to the J function, for example J\v( r ) is identically 1 for Poisson processes. This enables 
direct interpretation of uncorrected estimates of J, something not possible with uncor- 
rected estimates of either F, G or K. We propose a Monte Carlo test for complete spatial 
randomness based on testing whether Jw( r ) = 1- Computer simulations suggest this test 
5_i ' is at least as powerful as tests based on edge corrected estimators of J. 
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1 Introduction 

A spatial point pattern is often studied by estimating point process characteristics such as 
the empty space function F, the nearest-neighbour distance distribution function G, and the 
i^-function. Here we consider 

as advocated by |van Lieshout and Baddeley (19"9lS| ). The J function is identically equal 
to 1 for a Poisson process, and values of J(r) less than or greater than 1 are suggestive of 
clustering or regularity, respectively. (Note that it is of course possible to find non-Poisson 
point processes for which J{r) = 1, as Pedford and van den Berg (19971) nave shown). 

In practice, observation of the point process is usually restricted to some bounded window 
W . As a consequence, estimation of the summary functions, which is based on the measurement 
of various distances of the point process, is hampered by the "edge effects" (bias and censoring) 
introduced by restricting observation of these distances to W . In order to counter edge effects 
it is necessary to apply some form of edge correction to the empirical estimates of the summary 
functions. For further details see |Baddeley (1999| ); |Stoyan et al. (199% |Cressie (1991| ); 



Hipley ( 1988 



Unbiasedness is highly desirable when a summary function estimate is to be compared 
directly to the corresponding theoretical value for a point process model. However, as Diggle 
argues in the discussion of |Ripley (1977|) and in plGGLE (1983|) , unbiasedness is not essential 



when using a summary function estimator as the test statistic in a hypothesis test, since the 
bias will be accounted for in the null distribution of the test statistic. 

This paper studies the uncorrected estimator of J obtained by ignoring edge effects and 
computing J = (1 — G)/(l — F) from the uncorrected, empirical distributions G and F of 
distances observed in a compact window. It was prompted by the accidental discovery that 
this uncorrected estimator remarkably still yields values J(r) approximately equal to 1 for 
the Poisson process. An intuitive explanation is that the relative bias due to edge effects is 
roughly equal for the estimates of 1 — G and 1 — F, so that these biases approximately cancel 
in the ratio estimator of J. It follows that the uncorrected estimate of J could be used for the 
direct visual assessment of deviations from the Poisson process, something not possible with 
the uncorrected estimates of F, G or K. Our aim is to formalise this uncorrected estimator 
of J, and to investigate its use as a summary statistic and as a test statistic in point pattern 
analysis. 

The paper is organised as follows. Section 2 outlines and generalises some fundamental 
point process ideas. In Section 3 we define the Jw function estimated by the uncorrected 
procedure described above and derives some of its properties. In Section 4 we verify that the 
natural estimator of Jw is the uncorrected estimate Jw- Finally, Section 5 presents the results 
of a computational experiment to compare the power of Monte Carlo Tests constructed from 
estimates of J and Jw as well as estimation results from the simulations of various point process 
models in square and rectangular windows in R 2 and in a cubic window in R 3 . 

2 Background 

Let A be a stationary point process in H d with intensity A. (For details of the theory of point 
processes see |Daley and Vere-Jones, 1988fc |Cressie, 199 It [Stoyan et al., 1995J) . The 



empty space function F(r) is the probability of finding a point of the process within a radius 
r of an arbitrary fixed point: 

F(r) = P(In5(0,r)/0), (2) 

where B(x,r) denotes the ball of radius r centred at x. The nearest neighbour distance dis- 
tribution function G(r) is the probability of finding another point of the process in the ball of 
radius r centred at a "typical" point of the process: 

G(r)=P o! (ln5(0,r)^), (3) 

where P 0! denotes the reduced Palm distribution at the origin 0. Roughly speaking P 0! is the 
distribution of rest of the process X \ {0} given there is a point of the process at the origin (see 



Daley and Vere-Jones, 1988 



Let W be a compact observation window in H d with nonempty interior. The construction 
of estimators for F is based on the stationarity of X yielding 

F(r) = -^— [p(XnB(x,r)^<&)dx, (4) 



IV 



where \W\ denotes the volume of W. For G the Campbell- Mecke formula ( JStoyan et al., 1995| . 
(4.4.3)) gives 

G(r) = -i-E Y, HXnB(x,r)\{x}^®}. (5) 

' ' xexnw 

In both cases we need information about X fl B(x,r) for x G W, whereas we only observe 
X fl B(x,r) fl W. Usually this edge effect problem is countered by restricting the integration 
in (Q) and summation in (|J) to those points x for which B(x,r) C W (the "border method") 
or by weighting the contributions to the integral and sum so as to correct for the bias (see for 
example Paddeley, 1999|) . 

The uncorrected, empirical distributions of distances observed in the window W correspond 
to simply replacing X by X fl W in ([|) and (|5|). In order to investigate the effect of this, we 
extend F and G to functionals as follows. 

Definition 1 . For every compact set K C H d containing the origin define 

F(K) := P(XnK^([>) 
<&(K) := P 0! (XH K^0>) 
1 - <G(K) 



3(K) 



¥(K)- 



(Note that the empty space functional bears some relation to the contact distribution function 
QStoyan et al., 19951 , P- 105) in that H B {r) = F(rB)). ^From these we are able to define the 
window based J function. 



3 The Jw function 

Definition 2. For every compact set W C R d with nonempty interior let 

f[l-G(B(0,r)nW- x )]dx 

Jw{r) := J[l-F(B(0,r)nW- x )]dx (6) 

w 

be the window based J function, where W_ x = {y — x : y e PF} is the translate of W by 

-x e R d . 



If X is a stationary Poisson process, then by Slivnyak's Theorem ( [Stoyan et al., 1995| . 



(4.4.7)) P = P 0! . Thus F = G and we arrive at the following proposition. 
Proposition 1. Let X be a stationary Poisson process. Then 

Jw{ r ) = 1 f or oil W and r > 0. 

Explicit evaluation of Jw for other point process models seems difficult. However, we can 
show that Jw behaves similarly to the J function for ordered and clustered processes, suggesting 
it can also be interpreted in the same way as the J function. For some processes it is also possible 
to demonstrate that the Jw function exhibits less deviation from the Poisson hypothesis than 
the equivalent J function. 

Proposition 2. Suppose X is a process which is "ordered" in the sense that its J functional 
is non- decreasing, that is K\ C i^ 2 implies J{K\) < J (^2)- Then 

1 < Jwij) < Jij') f or oil r. 

Similarly if a process is "clustered" , K\ C K 2 implies J(Ki) > J(K 2 ), then 1 > Jw( r ) > J{ r ) 
for all r. 

PROOF: Observe that @ can be rewritten 

Jw(r)= I J{B(0,r)r\W- x )h w , r (x)dx (7) 

Jw 

where 

(l-F(B(0,r)nW. x )) 

WA) f w (l-HB(0,r)nW- y ))dy 
satisfies h w ,r(%) > for all x E H d and J w hw,r(x)dx = 1. Hence 

minJ(S(0,r)nW_ s ) < J w {r) < max J(B(0,r) n W- x ) 

w w 

and since J is nondecreasing 

1 = J(0) < J w ( r ) < J(5(0,r)) = J(r). D 
The latter result can be strengthened to strict inequality for specific examples. 



Proposition 3. Let X be a Neyman- Scott cluster process with mean number of points per 
cluster greater than 1. Assuming the support of the distribution of the cluster points contains 
a neighbourhood of the origin, then 

J(r) < J w (r) < 1 for all W and r > 0. 

Examples of processes satisfying the conditions of Proposition |3] are Matern's cluster process 
and the modified Thomas process described, for example, in [Stoyan et al. (1995[ ). 



PROOF: The Palm distribution of a Neyman-Scott process is the convolution of the original 
distribution P of the process and the Palm distribution cq of the representative cluster N 
(|Stoyan et al., 1995| , (5.3.2)). Thus, for every compact set K, we have 




l-G(K) = / l{(ip\Jxp)nK\{0} = ®}co(dip)P(d<p) 




Hence 



t{(p n K \ {0} = 0}i{^ n K \ {0} = 0}co(#)P(dp) 

p(x n k \ {0} = 0)c o (iv n k \ {0} = 0) 

(l-F(K))co(NnK\{0} = Q)). 



Now the assumption on the cluster distribution ensures cq(N fl K \ {0} = 0) < 1 for all K 
containing a neighbourhood of the origin. The conclusion then follows since -6(0, r) D W_ x 
contains a neighbourhood of whenever x is an interior point of W. □ 

Proposition 4. Let X be a hard-core process with hard-core radius R. Then 

1 < Jw(r) < J{r) for all W and < r < R, 
and Jw{ r ) is non- decreasing in r for all r < R. 

PROOF: Trivially we have G(K) = for all K C B(0,R), while for any point process the 
empty space functional F is nondecreasing. Therefore the J functional becomes 

JW = 1 _ I {K) for all K<zB(0,R), 

which is also non-decreasing and the result follows by proposition ^. □ 

4 Estimation of the J w function 

Analogously to J we want to estimate Jw by the ratio of two estimators for the denominator 
and numerator in Definition 0. The stationarity of X and Fubini's Theorem yield 

J F(5(0, r) n W. x ) dx = E 1 1{X n B(x, r) n W + 0} dx, 

w w 



and so the denominator of (fj|) becomes 

Ml - F(B(0, r) n W- x )\ dx=\W\ 



1- 



liy 



■E|^npn^)©B(o,r) 



(8) 



w 



where © denotes Minkowski addition. 

Applying the Campbell-Mecke formula QStoyan et al., 199"5| , (4.4.3)) we find 



[ <S(B(0,r)nW- x )dx = fp o! (lnB(0,r)nf.^(J)(li 
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w 



\e J2 !{XnB(x,r)\{x}nW ^(b}. 



xexnw 



Let d(x, A) denote the Euclidean distance from a point x G H d to a set A C R d . The numerator 
of (Bl) can then be expressed as 



f[l-G{B(0,r)r\W- a) )]dx= \W\ 1 - —J— E V t{d{x,X n VT \ {x}) < r} 
The two results (B) and (0) allow uncorrected estimation of Jw(r) by 



(9) 



1- 



i 



#(XrW) 



X) lKx,Inlf\{x})<r} 



■AvO") := 



xexnlf 



i-^|vyn((xnw^)©5(o,r) 



(10) 



which is the uncorrected estimate of the J function referred to in the introduction. 

Thus the uncorrected estimate of the J function, based on the uncorrected (EDF) estimates 
of F and G, can be thought of as a ratio unbiased estimator of the Jw function. As was shown 
in the previous section, the Jw function can be interpreted in the same way as the J function. 
Consequently, the uncorrected estimate of the J function, unlike the uncorrected estimates of 
F, G or K, can be used directly as an interpretive statistic in classifying deviations from the 
Poisson process. 



5 Simulation study 

This section reports the results of simulation studies comparing the uncorrected estimator Jw 



with "corrected" estimators J. In § |5.1| we show the results for a single simulated pattern; § p72 
reports the means and variances of the estimators of J in a simulation study. These results 
show that Jw typically has smaller variance than the corrected estimators. In §^]3| and § [5.4| we 
consider the power of hypothesis tests based on the uncorrected estimator J w . It is not clear, a 
priori, whether "corrected" estimators J or uncorrected estimators Jw will yield more powerful 
tests. There are two competing effects: the variance of Jw is smaller than that of J, but J w is 
less sensitive than J to departures from the Poisson process (say) according to Proposition |2|. 
For the comparisons which follow, the reduced sample (border method) estimator J rs was 
adopted. However, for the purpose of highlighting the variations which exist between corrected 



estimates of the J function, the Kaplan-Meier estimator Jk m (described in Paddeley and 
Pill, 1997] ) is also considered in many cases. 



5.1 Empirical example 

This example highlights the use of Jw as a qualitative summary statistic. Consider the point 
pattern given in Figure |l|. This pattern is a realisation of a Matern Cluster Process, intensity 
A = 100, cluster radius R = 0.1 and mean number of offspring n = 4, observed within an 
observation window consisting of two rectangular windows, 3.125 by 0.16 units, separated 
by 0.02 units. (The Matern Cluster Process is a Neyman-Scott process in which offspring 
are uniformly distributed in disc of radius R about (Poisson) parent points. The number of 
offspring per parent point is Poisson with mean fi QStoyan et al., 1995| , p. 159)). 

Figure also displays the corresponding estimates Jw and J rs of the given point pattern, 
together with envelopes of 99 simulations of a binomial process of the same intensity, observed 
within the same window. The results in this case are marked; the uncorrected estimate suggests 
strong evidence of clustering in the point pattern, while the corrected estimate appears to 
suggest no evidence of clustering. Of course, the results are not surprising given the severity of 
the edge effect introduced by the "censoring" of the middle seventeenth of the window compared 
to the relatively small bias this introduces. However, it does illustrate the possible benefit of 
using Jw in certain situations. 

As the empirical use of Jw is the same as the J function, readers interested in further 
examples of the analysis of empirical data are referred to Kerscher (1998 ), Kerscher et al.' 
KT999D , and |Kerscher et al. (1998|) . 



5.2 Mean and variance 

^From the previous example it is clear that the uncorrected estimate of the J function may be 
superior in some situations. To examine whether this was true more generally, a number of 
simulations were conducted to compare the corrected and uncorrected methods across a range 
of processes. This began with the estimation of the mean and standard deviation of the three 
estimators ( Jw, Jkm, and J rs ) based on 10, 000 realisations of a Poisson process with intensity 
100, in a unit square window, with the results presented in Figure |2|. 

With increasing r, the distributions of the J (that is, the estimators Jw, Jkrn, and J rs ) 
become skewed to the right; for large r there is substantial mass above J.( r ) = 2. As a 
result all three estimators are positively biased for large values of r. Empirically it was found 
that a square root transformation approximately symmetrised the distribution. As expected, 
the sample standard deviation of the estimates increases with r, as the denominator of each 
estimator decreases with r. However, Jw is less biased and has lower variance than Jkm and 
j rs . 

These simulations were repeated for two processes with more substantial edge effects (namely, 
a Poisson process of intensity A = 25 in a unit window, and a Poisson process of intensity A = 10 
in a 10 by 1 rectangular window). In both cases the results were qualitatively similar to those 
above. 

In addition, further simulations were conducted for point patterns in R 3 . Estimates of the 
means and standard deviations of Jk m and J rs based on 1000 realisations in a unit cube were 



compared for the Poisson process and two alternatives: Matern hard-core QStoyan et al., 1995| , 
p. 163) and Matern cluster processes for a range of parameter values. Some of the key results 
are presented in Figure [| As expected, Jw is reliable over a wider domain than J rs . For 
hard-core processes, the standard deviation of J rs is considerably bigger than that of Jw and 



the difference grows with the hard-core radius. For cluster processes the differences are far less 
apparent, however the overall tendency of Jw to have lower variance is also confirmed for this 
class of processes. Note also that with both alternative processes the mean of Jw is bounded 
by 1 and J rs , as expected. 

It is interesting to note that, unlike the J function estimators, the domain of the Jw function 
estimator for a given point process realisation can be easily calculated. The Jw estimator is 
defined for all r < rp max , where rp max is the maximum nearest-point distance (the largest 
distance, over all points in the window, from a point to the nearest point of the process). Also 
Jw(t) = for any rn < r < rp if rn < rp , where ra is the maximum nearest- 

yy V / J ^max — J^max ^Jmax ? max * ^max 

neighbour distance (the largest distance, over all points of the process within the window, from 
a point of the process to the nearest other point of the process). The value rp max is however an 
upper bound on the domain of both the Reduced Sample and Kaplan-Meier estimators. 

5.3 The test statistic 

We now aim to compare the power of the Jw function with the edge corrected estimators of 
the J function in testing the Poisson hypothesis in the two dimensional case. We restricted 
ourselves to this estimation problem in view of the problems with estimating the range of 
interaction using the J function reported in [Kerscher et al. (1999| ) . 



The distribution of the following test statistic for each of the three estimators was estimated: 

" r ° J(r) - 1 



■dr, (11) 



a{r) 

where a denotes the sample standard deviation of J (r) under the Poisson hypothesis. This form 
of test statistic was chosen, as opposed to a squared integrand, because of the skewed nature 
of the distributions of J. The distributions of the test statistics were estimated by a discrete 
sum and based on 10, 000 realisations of a Poisson process. The upper limit of integration r 
was chosen to be the 0.9 quantile of the F function (for intensity 100, r$ p» 0.856). Having 
estimated the distribution, the 0.025 and 0.975 quantiles were obtained for use in a two-sided 
5% significance test for deviation from a Poisson process. One-sided 5% significance tests were 
also constructed to test for clustering or regularity by considering the 0.05 and 0.95 quantiles, 
respectively. 

5.4 Power of tests using the various estimators 

In order to estimate the power of the hypothesis test described above, realisations from alter- 
native point processes were generated and the proportion of the realisations rejected by the 
hypothesis test was recorded. The first class of point processes considered was the Matern 
hard-core process, with hard-core radius R. For each of 22 values of R, 1000 realisations were 
generated. The proportion of rejections is presented in Figure |. Note that as R — > the model 
approaches the Poisson process, so we expect all power curves to approach 0.05 (5%) as R — > 0. 
All three estimators have very similar power curves, with the Jw estimator at least as powerful 
as the two J function estimators for all values of R. 

The other class of alternative point processes considered was Matern's cluster process. A 
grid of (R, n) values was constructed and 1000 realisations were obtained from the corresponding 
Matern cluster process. The proportion of rejections for each (R, //) value is presented in Figure 

8 



g| Once again, the curves are very similar, with all three tests performing almost identically. 
Similar results were obtained for the respective one-sided tests in both cases as well as for the 
lower intensity 25. 

The power tests for the 10 by 1 window support the argument that edge effects are stronger 
when the boundary is relatively longer. The resulting power function estimates against the 
Matern Model II and Matern cluster process models are also presented in Figures |] and |^, 
respectively. 

One important observation made while conducting these numerical simulations was that 
the choice of test statistic had far more impact on the power of the resulting hypothesis test 
than the choice of J function estimator. For a comparison of various test statistics of the J 
function see |Thonnes and van Lieshout (1998Q . 
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Figure 1: Top: empirical data. Bottom: empirical J rs (left) and Jw (right) functions (points) 
and envelope of 99 simulations of a binomial process with the same intensity (solid lines). 
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Figure 2: Results in the unit square, intensity 100. Mean (top) and standard deviation (bottom) 
of J estimators as a function of r for a Poisson process 
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Figure 3: Results in the unit cube, intensity=100. Top: Mean (left) and standard deviation 
(right) of J estimators as a function of r for a Poisson process. Bottom left: Mean of J 
estimators for a Matern model II process (R = 0.1). Bottom right: Mean of J estimators for a 
Matern cluster process (/i — 4, R — 0.1). 
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Figure 4: Power against Matern Model II as a function of hard-core radius R. Top: unit square, 
intensity 100. Bottom: [0, 1] x [0, 10] rectangle, intensity 10 
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Figure 5: Power against Matern Cluster Process as a function of mean cluster size for various 
cluster radii R. Top: unit square, intensity 100. Bottom: [0, 1] x [0, 10] rectangle, intensity 10 
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